﻿using System;
using System.Text;
using System.Drawing;
using System.Buffers;
using System.Collections;
using System.Collections.Generic;
using System.Runtime.InteropServices;

public static partial class NativeAOT
{
    [UnmanagedCallersOnly(EntryPoint = "remz")]
    public static unsafe double remz(double a, double b, IntPtr p_ptr, int n, double eps, IntPtr f_x_ptr)
    {
        double* p = (double*)p_ptr.ToPointer();
        f_x = Marshal.GetDelegateForFunctionPointer<delegatefunc_x>(f_x_ptr);

        return remz(a, b, p, n, eps);
    }

    /// <summary>
    /// Remez算法
    /// f计算f(x)函数值的函数名。
    /// </summary>
    /// <param name="a">区间左端点值。</param>
    /// <param name="b">区间右端点值。</param>
    /// <param name="p">p[n]返回n-1次最佳一致逼近多项式的系数。</param>
    /// <param name="n">n-1次最佳一致逼近多项式的项数。</param>
    /// <param name="eps">控制精度要求。</param>
    /// <returns>函数返回偏差绝对值。</returns>
    public static unsafe double remz(double a, double b, double* p, int n, double eps)
    {
        int i, j, k, m;
        double d, t, u, s, xx, x0, h, yy;
        double[] x = new double[21];
        double[] g = new double[21];
        double[] pp = new double[21];

        if (n > 20) n = 20;
        m = n + 1;
        d = 1.0e+35;

        // 初始点集
        for (k = 0; k <= n; k++)
        {
            t = Math.Cos((n - k) * 3.1415926 / (1.0 * n));
            x[k] = (b + a + (b - a) * t) / 2.0;
        }

        while (true)
        {
            u = 1.0;
            for (i = 0; i <= m - 1; i++)
            {
                pp[i] = f_x(x[i]);
                g[i] = -u; u = -u;
            }
            for (j = 0; j <= n - 1; j++)
            {
                k = m;
                s = pp[k - 1];
                xx = g[k - 1];
                for (i = j; i <= n - 1; i++)
                {
                    t = pp[n - i + j - 1];
                    x0 = g[n - i + j - 1];
                    pp[k - 1] = (s - t) / (x[k - 1] - x[m - i - 2]);
                    g[k - 1] = (xx - x0) / (x[k - 1] - x[m - i - 2]);
                    k = n - i + j;
                    s = t;
                    xx = x0;
                }
            }
            u = -pp[m - 1] / g[m - 1];
            for (i = 0; i <= m - 1; i++)
            {
                pp[i] = pp[i] + g[i] * u;
            }
            for (j = 1; j <= n - 1; j++)
            {
                k = n - j;
                h = x[k - 1];
                s = pp[k - 1];
                for (i = m - j; i <= n; i++)
                {
                    t = pp[i - 1];
                    pp[k - 1] = s - h * t;
                    s = t;
                    k = i;
                }
            }
            pp[m - 1] = Math.Abs(u);
            u = pp[m - 1];
            if (Math.Abs(u - d) <= eps)
            {
                for (i = 0; i < n; i++)
                {
                    p[i] = pp[i];
                }
                return (u);
            }
            d = u;
            h = 0.1 * (b - a) / (1.0 * n);
            xx = a;
            x0 = a;
            while (x0 <= b)
            {
                s = f_x(x0);
                t = pp[n - 1];
                for (i = n - 2; i >= 0; i--)
                {
                    t = t * x0 + pp[i];
                }
                s = Math.Abs(s - t);
                if (s > u)
                {
                    u = s;
                    xx = x0;
                }
                x0 = x0 + h;
            }
            s = f_x(xx);
            t = pp[n - 1];
            for (i = n - 2; i >= 0; i--)
            {
                t = t * xx + pp[i];
            }
            yy = s - t;
            i = 1;
            j = n + 1;
            while ((j - i) != 1)
            {
                k = (i + j) / 2;
                if (xx < x[k - 1]) j = k;
                else i = k;
            }
            if (xx < x[0])
            {
                s = f_x(x[0]);
                t = pp[n - 1];
                for (k = n - 2; k >= 0; k--)
                {
                    t = t * x[0] + pp[k];
                }
                s = s - t;
                if (s * yy > 0.0)
                {
                    x[0] = xx;
                }
                else
                {
                    for (k = n - 1; k >= 0; k--)
                    {
                        x[k + 1] = x[k];
                    }
                    x[0] = xx;
                }
            }
            else
            {
                if (xx > x[n])
                {
                    s = f_x(x[n]);
                    t = pp[n - 1];
                    for (k = n - 2; k >= 0; k--)
                    {
                        t = t * x[n] + pp[k];
                    }
                    s = s - t;
                    if (s * yy > 0.0)
                    {
                        x[n] = xx;
                    }
                    else
                    {
                        for (k = 0; k <= n - 1; k++)
                        {
                            x[k] = x[k + 1];
                        }
                        x[n] = xx;
                    }
                }
                else
                {
                    i = i - 1;
                    j = j - 1;
                    s = f_x(x[i]);
                    t = pp[n - 1];
                    for (k = n - 2; k >= 0; k--)
                    {
                        t = t * x[i] + pp[k];
                    }
                    s = s - t;
                    if (s * yy > 0.0) x[i] = xx;
                    else x[j] = xx;
                }
            }
        }
    }

    /*
    // Remez算法例
      int main()
      { 
          int i;
          double a,b,eps,p[4], u;
          double remzf(double);
          a=-1.0; b=1.0; eps=1.0e-10;
          u = remz(a,b,p,4,eps,remzf);
          cout <<"最佳一致逼近多项式系数 :" <<endl;
          for (i=0; i<=3; i++)
             cout <<"p(" <<i <<") = " <<p[i] <<endl;
          cout <<"偏差绝对值 = " <<u <<endl;
          return 0;
      }

      double remzf(double x)
      { 
          return(exp(x));
      }
    */
}

